DATA 605 - Final Project

Exploring the quantitative variables

kaggle_numerics %>% 
  dplyr::select(2:11) %>%
  pairs.panels(method = "pearson", hist.col = "#8bb397")

kaggle_numerics %>% 
  dplyr::select(12:21) %>%
  pairs.panels(method = "pearson", hist.col = "#8bb397")

kaggle_numerics %>% 
  dplyr::select(22:31) %>%
  pairs.panels(method = "pearson", hist.col = "#8bb397")

kaggle_numerics %>% 
  dplyr::select(32:38) %>%
  pairs.panels(method = "pearson", hist.col = "#8bb397")

After creating a few plot matrices with our quantitative data, I thought there were a few initial takeaways:

• Many of the living spaces, and the square footage values seem to show stronger correlations with one another. This makes sense given that many times, the square footage of one floor will likely be similar square footage to another floor in the house. For instance, there’s a strong correlation between the first floor square footage and the basement square footage. When thinking about our regression later, it’s nice to know that there are multiple attributes here that could serve as proxies (if needed).

• Additionally, I’m seeing a gradual, positive trend towards larger areas of garages the later they are built in the 20th century – for instance, garages built in the early part of the 20th century seem to be smaller than those that are built in the later 1900s and early 2000s. This may be something to look into further, and it’s relationship with the Sale Price later on.

• It’s interesting to see too that the year the house is built has a relationship with the rating of overall material and finish quality of the house. Again, this make sense, but it would be valueable to see if this has a relationship with Sale Price.

Now, let’s explore the categorical variables a bit.

Exploring a few of the categorical variables

table(kaggle_categorical$LandSlope)

## 
##  Gtl  Mod  Sev 
## 1382   65   13

# table(kaggle_categorical$BldgType)
table(kaggle_categorical$HouseStyle)

## 
## 1.5Fin 1.5Unf 1Story 2.5Fin 2.5Unf 2Story SFoyer   SLvl 
##    154     14    726      8     11    445     37     65

# table(kaggle_categorical$RoofStyle)
table(kaggle_categorical$RoofMatl)

## 
## ClyTile CompShg Membran   Metal    Roll Tar&Grv WdShake WdShngl 
##       1    1434       1       1       1      11       5       6

# table(kaggle_categorical$ExterCond)
table(kaggle_categorical$CentralAir)

## 
##    N    Y 
##   95 1365

# table(kaggle_categorical$Utilities)
table(kaggle_categorical$HeatingQC)

## 
##  Ex  Fa  Gd  Po  TA 
## 741  49 241   1 428

table(kaggle_categorical$SaleType)

## 
##   COD   Con ConLD ConLI ConLw   CWD   New   Oth    WD 
##    43     2     9     5     5     4   122     3  1267

table(kaggle_categorical$SaleCondition)

## 
## Abnorml AdjLand  Alloca  Family  Normal Partial 
##     101       4      12      20    1198     125

table(kaggle_categorical$Foundation)

## 
## BrkTil CBlock  PConc   Slab  Stone   Wood 
##    146    634    647     24      6      3

After looking through a large number of categorical values, and checking frequencies, a few things stood out to me:

• It looks like there’s a majority of houses that are built on gentle slopes, however, there’s a proportion that are built on moderate or severe slopes – this could be something to look into further when building out the regression later on. It may be interesting to see if there is any correlation here with the Sale Price.

• Similarly, the foundation material seems to be split heavily between cinder blocks and poured concrete. This is another feature that may be useful to explore later to see if this has an affect on housing prices, or if there is a connection between the year a house was built and which type of foundation was used.

• Finally, Sale Type and Sale Condition were interesting values and their frequencies suggest that there isn’t a “one-size-fits-all” approach to certain transactions, making this an enticing set of features to explore further.

Developing a scatterplot matrix

After plotting all of my quantitative variables earlier, I thought I’d hone in on a few that may be useful to explore their relationship with the Sale Price variable. You can find my plot below:

pairs(kaggle_df[, c('1stFlrSF', 'GrLivArea', 'GarageArea', 'TotalBsmtSF', 'YearBuilt', 'SalePrice' )], pch = 19)

Developing a correlation matrix

For our correlation matrix, I will use GrLivArea, 1stFlrSF, and SalePrice. And taking a few of these same attributes, I thought it would be interesting to plot these in a correlation matrix:

kaggle_numerics %>% 
  dplyr::select(`1stFlrSF`, `GrLivArea`, SalePrice) %>%
  pairs.panels(method = "pearson", hist.col = "#8bb397")

corr_matrix <- kaggle_numerics %>% 
  select(`1stFlrSF`, GrLivArea, SalePrice) %>% 
  cor() %>% 
  as.matrix()

As we can see from above, all three attributes, the square footage of the first floor, the above grade (ground) living area and the Sale Price all seem to have a positive relationship with one another. All three distributions also appear to be unimodal and right-skewed. Below is the official correlation matrix:

kable(corr_matrix) %>% 
  kable_styling(bootstrap_options = 'striped', full_width = FALSE)

Hypothesis testing on correlation matrix

Now, although these correlations seem to be quite strong, we can do some hypothesis testing on their relationships at an 80-percent confidence interval:

corr1 <- cor.test(kaggle_numerics$SalePrice, kaggle_numerics$`1stFlrSF`, method = 'pearson', conf.level = 0.80)

corr2 <- cor.test(kaggle_numerics$SalePrice, kaggle_numerics$GrLivArea, method = 'pearson', conf.level = 0.80)

corr3 <- cor.test(kaggle_numerics$`1stFlrSF`, kaggle_numerics$GrLivArea, method = 'pearson', conf.level = 0.80)

corr1

## 
##  Pearson's product-moment correlation
## 
## data:  kaggle_numerics$SalePrice and kaggle_numerics$`1stFlrSF`
## t = 29.078, df = 1458, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 80 percent confidence interval:
##  0.5841687 0.6266715
## sample estimates:
##       cor 
## 0.6058522

corr2

## 
##  Pearson's product-moment correlation
## 
## data:  kaggle_numerics$SalePrice and kaggle_numerics$GrLivArea
## t = 38.348, df = 1458, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 80 percent confidence interval:
##  0.6915087 0.7249450
## sample estimates:
##       cor 
## 0.7086245

corr3

## 
##  Pearson's product-moment correlation
## 
## data:  kaggle_numerics$`1stFlrSF` and kaggle_numerics$GrLivArea
## t = 26.217, df = 1458, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 80 percent confidence interval:
##  0.5427732 0.5884078
## sample estimates:
##      cor 
## 0.566024

Solution: After running the three pairwise correlations between my variables of SalePrice, 1stFlrSF, and GrLivArea, I can see from the outputs that all three accept the alternative hypothesis that the true correlation is not equal to zero and that there is a positive relationship between these attributes. Therefore, this rejects the null hypothesis that states that the correlations between each pairwise set of variables is zero. I know this to be true since the p-values for all three of my Person’s correlation tests are less than 0.05.

Solution: Additionally, since we are dealing with a large sample size, and our correlation tests are yielding very small p-values, all less than 0.001, I’m not worried about family-wise error when measuring relationships across these three attributes. With p-values this low, it creates a solid argument that we aren’t running into a Type 1 error here.

Linear Algebra and Correlation: Invert your correlation matrix from above. (This is known as the precision matrix and contains variance inflation factors on the diagonal.) Multiply the correlation matrix by the precision matrix, and then multiply the precision matrix by the correlation matrix. Conduct LU decomposition on the matrix.

Remembering that the inverse of a matrix multiplied by a matrix will become the identity matrix:

\[A^{−1}A=AA^{−1}=I\]

We can find our precision matrix by using the solve() function in baseR to invert our correlation matrix:

prec_matrix <- solve(corr_matrix)
prec_matrix

##             1stFlrSF  GrLivArea SalePrice
## 1stFlrSF   1.6795240 -0.4611713 -0.690746
## GrLivArea -0.4611713  2.1352622 -1.233697
## SalePrice -0.6907460 -1.2336974  2.292718

Then, we can multiple the precision matrix by the correlation matrix to confirm that we get the identity matrix:

\[ correlation \ matrix \times \ precision \ matrix = \]

I_1 <- round(prec_matrix %*% corr_matrix, 2)
I_2 <- round(corr_matrix %*% prec_matrix, 2)

I_1

##           1stFlrSF GrLivArea SalePrice
## 1stFlrSF         1         0         0
## GrLivArea        0         1         0
## SalePrice        0         0         1

I_2

##           1stFlrSF GrLivArea SalePrice
## 1stFlrSF         1         0         0
## GrLivArea        0         1         0
## SalePrice        0         0         1

As we can see above, we do indeed get the Identity Matrix when we multiply the precision matrix with the correlation matrix and vice-versa.

Finally, we will perform LU decomposition on the matrix. Fortunately, R has a package called lu.decomposition() that will find the upper and lower triangle matrices of our precision matrix, \(U\) and \(L\) respectively. The relationship between these two matrices is based on the formula:

\[A = LU\] Where \(A\) = our precision matrix, and \(L\) and \(U\) are the upper and lower triangular matrices of our precision matrix. We can see, when finding \(L\) and \(U\) below, we get these respective matrices:

A <- prec_matrix
L <- lu.decomposition(A)$L
U <- lu.decomposition(A)$U

round(L, 2)

##       [,1]  [,2] [,3]
## [1,]  1.00  0.00    0
## [2,] -0.27  1.00    0
## [3,] -0.41 -0.71    1

round(U, 2)

##      [,1]  [,2]  [,3]
## [1,] 1.68 -0.46 -0.69
## [2,] 0.00  2.01 -1.42
## [3,] 0.00  0.00  1.00

And when we multiply these together, we can obtain our original precision matrix (\(A\)):

round(A, 2)

##           1stFlrSF GrLivArea SalePrice
## 1stFlrSF      1.68     -0.46     -0.69
## GrLivArea    -0.46      2.14     -1.23
## SalePrice    -0.69     -1.23      2.29

round(L %*% U, 2)

##       [,1]  [,2]  [,3]
## [1,]  1.68 -0.46 -0.69
## [2,] -0.46  2.14 -1.23
## [3,] -0.69 -1.23  2.29

\[\begin{equation} LU = A, where \space A = \begin{bmatrix} 1.68 & -0.46 & -0.69\\ -0.46 & 2.14 & -1.23\\ -0.69 & -1.23 & 2.29\\ \end{bmatrix} \space = \space \begin{bmatrix} 1.00 & 0.00 & 0.00\\ -0.27 & 1.00 & 0.00\\ -0.41 & -0.71 & 1.00\\ \end{bmatrix} \space \begin{bmatrix} 1.68 & -0.46 & -0.69\\ 0.00 & 2.01 & -1.42\\ 0.00 & 0.00 & 1.00\\ \end{bmatrix} \space = LU \end{equation}\]

Calculus-Based Probability & Statistics: Many times, it makes sense to fit a closed form distribution to data. Select a variable in the Kaggle.com training dataset that is skewed to the right, shift it so that the minimum value is absolutely above zero if necessary. Then load the MASS package and run fitdistr to fit an exponential probability density function. (See https://stat.ethz.ch/R-manual/R-devel/library/MASS/html/fitdistr.html ). Find the optimal value of \(\lambda\) for this distribution, and then take 1000 samples from this exponential distribution using this value (e.g., rexp(1000, \(\lambda\))). Plot a histogram and compare it with a histogram of your original variable. Using the exponential pdf, find the 5th and 95th percentiles using the cumulative distribution function (CDF). Also generate a 95% confidence interval from the empirical data, assuming normality. Finally, provide the empirical 5th percentile and 95th percentile of the data. Discuss.

For this exercise, we’ll use the GrLivArea variable from the training dataset. When we plot it below on a histogram we can quickly see that it is skewed to the right:

ggplot(kaggle_df, aes(GrLivArea)) +
  geom_histogram(aes(y = ..density..), colour="#555555", fill="#dddddd") +
  geom_density(alpha=.6, fill="#333333") +
  geom_vline(xintercept = min(kaggle_df$GrLivArea)) +
  labs(title = "Above Ground Living Area (square feet)", align='middle') +
  theme_minimal()

We can also see that the minimum value is absolutely above zero, so we do not need to perform any shifts. We can double check the data by using the summary() function, and we see that the minimum value is 334 square feet:

summary(kaggle_df$GrLivArea)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     334    1130    1464    1515    1777    5642

Loading the MASS package and running fitdistr()

After loading the MASS package and running fitdistr() to fit an exponential probability density function to our variable, I’ve found that that the optimal value of \(\lambda\) for this distribution is 0.000659864.

require(MASS)

exp_fit <- fitdistr(kaggle_df$GrLivArea, "exponential")
est_rate <- exp_fit$estimate
est_rate

##        rate 
## 0.000659864

Plot a histogram comparing the exponential distribution to the original

Using my newly found optimal value for lambda, I then was able to generate a 1000-sample variable from this exponential distribution:

set.seed(123)
pdf <- rexp(1000, est_rate)
pdf_df <- as.data.frame(pdf)

I then plotted this distribution below, along with my original distribution for GrLivArea:

A few takeaways when comparing the distributions:

• The original data for GrLivArea appears to be shaped as a more normal distribution, albeit right skewed, the most frequent square footage tends to be around 1500 and 1650 square feet.

• Alternatively, the exponential distribution of GrLiveArea appears to be even more right skewed, and takes on more of an exponential shape with a longer tail.

Find the 5th and 95th percentiles using CDF

Next, I’ll find the 5th and 95th percentiles using the cumulative distribution function (CDF). Thinking back to the formula from our reading, this function is:

\[F(x) = 1 - e^{-\lambda x}\] Since we already have found \(\lambda = 0.000659864\), we can solve the equation accordingly for the 5th percentile:

\[0.05 = 1 - e^{-0.000659864x}\] \[-0.95 = e^{-0.000659864x}\] \[-ln(0.95) = 0.000659864x\] \[\frac{-ln(0.95)}{0.000659864} = x\] \[77.73313 = x\] Similarly, we can find the 95th percentile by:

\[0.95 = 1 - e^{-0.000659864x}\] \[-0.05 = e^{-0.000659864x}\] \[-ln(0.05) = 0.000659864x\] \[\frac{-ln(0.05)}{0.000659864} = x\] \[4539.92 = x\] We can check these calculations using the qexp() function in R:

qexp(0.05, rate=est_rate)

## [1] 77.73313

qexp(0.95, rate=est_rate)

## [1] 4539.924

It looks like these match the calculations above!

Generate 95% confidence interval from empirical data

Assuming normality, we can generate a 95% confidence interval from the empirical data by utilizing the following function:

\[\overline{X} \pm Z * \frac{\sigma}{\sqrt{n}}\]

Since we have all of this information available to us, we can solve by utilizing the following:

• \(Z\) = 1.96 (standard for a 95% confidence interval)

• \(\overline{X}\) = mean value of GrLivArea

• \(\sigma\) = standard deviation of GrLivArea mean

• \(n\) = sample size of GrLivArea

X_cap <- mean(kaggle_df$GrLivArea)
Z <- 1.96
sigma <- std(kaggle_df$GrLivArea)
n <- length(kaggle_df$GrLivArea)

upper_bound <- X_cap + Z * (sigma / sqrt(n))
lower_bound <- X_cap - Z * (sigma / sqrt(n))

upper_bound

## [1] 1542.419

lower_bound

## [1] 1488.509

### talk about this tomorrow, notice that because original distribution is right skewed (not exactly normal), median would be better approximation here, so assuming normal distribution wouldn't be best here 1515.46 (mean) --> 1464 (median)

### CI [1542 - 1488], median falls outside of this boundary

mean(kaggle_df$GrLivArea)

## [1] 1515.464

median(kaggle_df$GrLivArea)

## [1] 1464

quantile(kaggle_df$GrLivArea, c(0.05, 0.95))

##     5%    95% 
##  848.0 2466.1